class: center, middle, inverse, title-slide .title[ # ADIM: Bayesian GLM. MCMC methods ] .subtitle[ ## Master’s Degree in Data Analysis, Process Improvement and Decision Support Engineering ] .author[ ###
Joaquín Martínez-Minaya, 2024-12-02
VAlencia BAyesian Research Group
Statistical Modeling Ecology Group
Grupo de Ingeniería Estadística Multivariante
] .date[ ###
jmarmin@eio.upv.es
] --- class: inverse, center, middle, animated, slideInRight # Motivation example --- # Heart Disease .left-column9[ - The study examines the relationship between: - the .hlb[myocardial infarction] (MI): `\(y = 1\)` if MI occurrence, or `\(y = 0\)` if No MI occurrence; and - .hlb[Age60]: Patients aged `\(\geq 60\)` (1) versus `\(< 60\)` (0). - .hlb[Systolic blood pressure (SBP140)]: SBP `\(\geq 140\)` mmHg (1) versus `\(< 140\)` mmHg (0). - .hlbred[Objective]: - Evaluate the association of age60 and sbp140 with MI probability. - Interpret the odds ratio (OR) for both predictors. ] -- .right-column9[ Table: Summary of Data from Study | y | age60 | sbp140 | |:-:|:-----:|:------:| | 0 | <60 | >=140 | | 0 | >=60 | <140 | | 0 | <60 | >=140 | | 0 | >=60 | >=140 | | 0 | >=60 | <140 | | 1 | <60 | <140 | ] --- # Bayesian Logistic Regression Model - .hlbpurple[Logistic regression] is used to model MI's probability based on age60 and sbp140. - .hlb[Likelihood] $$ y_i \sim \text{Bernoulli}(\pi_i)\,, i= 1, \ldots, 400\,, $$ using logit link: $$ logit(\pi_i) = \beta_0 + \beta_1 age60_i + \beta_2 sbp140_i$$ - .hlb[Prior distributions] (weakly-informative): \begin{eqnarray} \beta_0 & \sim & \mathcal{N}(0, 10^3), \nonumber \\ \beta_1 & \sim & \mathcal{N}(0, 10^3), \nonumber \\ \beta_2 & \sim & \mathcal{N}(0, 10^3), \nonumber \end{eqnarray} -- .hlbred[Note: There are no conjugate priors available for the logistic regression model.] --- #Table of contents ## 1. Bayesian computation. MCMC methods ## 2. Bayesian Software for MCMC --- class: center, middle, animated, rotateInUpRight, inverse # 1. Bayesian computation. MCMC methods --- # Monte Carlo Methods .left-column9[ .hlbpurple[Monte Carlo Simulation] - Draw .hlb[realizations of a random variable] for which only its density function is (fully or partially) known. ``` r x <- rnorm(1000, mean = 1, sd = 2) ``` <!-- --> ] -- .right-column9[ .hlbpurple[Monte Carlo Integration] - Computing the mean of a N(1, 2), - `\(E(X) = \int_{-\infty}^{\infty} x \cdot f(x) dx\)` - Using .hlb[Monte Carlo integration]: - Simulate from N(1, 2^2): `\(\phi^1, \ldots, \phi^N\)`. - Compute the mean of the simulated values: `\(E(X) \approx \frac{1}{N} \sum_{i=1}^N \phi^i\)` - Doing `summary` of the simulation, we compute more measures: | Min.| 1st Qu.| Median| Mean| 3rd Qu.| Max.| |------:|-------:|------:|----:|-------:|-----:| | -6.065| -0.334| 1| 1| 2.332| 8.086| ] --- # Markov Chain Monte Carlo - A Markov chain is a .hlbred[stochastic sequence of numbers] where each value in the sequence depends only upon the last. -- - If `\(\phi^1, \phi^2, \ldots, \phi^N\)` is a sequence of numbers, then `\(\phi^2\)` is only a function of `\(\phi^1\)`, `\(\phi^3\)` of `\(\phi^2\)`, etc. -- - Under certain conditions, the distribution over the states of the .hlb[Markov chain] will .hlb[converge to a stationary distribution]. -- - The .hlb[stationary distribution is independent of the initial starting values] specified for the chains. -- - AIM: construct a Markov chain such that .hlb[the stationary distribution is equal to the posterior distribution] `\(p(\theta \mid x)\)`. -- - We combine Markov Chain with Monte Carlo simulation --> .hlbpurple[Markov chain Monte Carlo (MCMC)]. -- - They were proposed by first time in the Statistics area by <a href = "https://www.tandfonline.com/doi/abs/10.1080/01621459.1990.10476213"> Gelfand and Smith (1990) </a>. --- # Posterior distribution ## Estimating the probability to score a penalty ### - .hlb[Likelihood] $$p(\boldsymbol{y} \mid \pi) = \pi^{k}(1-\pi)^{N-k} $$ -- ### - .hlb[Prior distribution] $$p(\pi) = \pi^{a-1}(1-\pi)^{b-1} $$ -- ### - .hlb[Posterior distribution] `$$p(\pi \mid \boldsymbol{y}) \propto p(\boldsymbol{y} \mid \pi) \times p(\pi) \propto \pi^{k + a - 1}(1 - \pi)^{N - k + b - 1} = p^*(\pi)$$` --- # MCMC: Metropolis-Hastings (MH) 1. Starting value `\(\pi^{(0)}\)` 2. For `\(t = 1, \ldots, T\)` - .hlb[We define a proposal distribution] (Usually similar to the objective distribution). In this case, `\(q(\pi \mid \pi^{(t-1)}) \sim logit-N(\pi^{(t-1)}, \sigma = 0.5)\)`. .hlb[Simulate] `\(\pi^{(prop)}\)` from it. - Compute .hlbred[probability of acceptance]: `$$\alpha = \text{min}\left(1, \frac{p^*(\pi^{(prop)})q(\pi^{(t-1)} \mid \pi^{(prop)})}{p^*(\pi^{(t-1)})q(\pi^{(prop)} \mid \pi^{(t-1)})}\right)$$` - Generate a .hlb[random number] `\(u\)` from the Uniform(0, 1). - `\(\pi^{(t+1)} = \pi^{(prop)}\)`, if `\(u \geq \alpha\)`, - `\(\pi^{(t+1)} = \pi^{(t)}\)`, if `\(u< \alpha\)` 3. Finally, we .hlb[obtain] `\(\pi^{0}, \pi^{1}, \ldots, \pi^{T}\)` which is .hlb[a simulation of the posterior distribution]. --- # MCMC: Metropolis-Hastings (MH) .center[  ] --- # Approaching probability of score using MH .left-column9[ .hlbpurple[Visual Metropolis-Hastings]
- Play the video from .hlb[minute 4:44]. ] -- .right-column9[ .hlb[Tracing the chain. Is the chain autocorrelated?]  ] --- # MCMC. Burnin and thin .left-column9[  ] .right-column9[   ] --- class: center, middle, animated, rotateInUpRight, inverse # 2. Bayesian Software for MCMC --- # Software .left-column9[ - <a href = "https://mcmc-jags.sourceforge.io/"> JAGS </a> - <a href="https://mc-stan.org/"> Stan </a> - <a href = "https://paulbuerkner.com/brms/"> brms </a>: allow for easy Bayesian Inference using Hamiltonian Monte Carlo - <a href="https://www.r-inla.org/"> INLA </a> (it does not use MCMC methods, but it is a very powerful tool) - <a href = "https://inlabru-org.github.io/inlabru/"> inlabru </a>: facilitate Bayesian Inference in Spatio-temporal models - <a href= "https://cran.r-project.org/web/packages/MCMCpack/index.html"> MCMC pack </a> (With this R-package, you can use MCMC methods with similar notation as usually use in R) - <a href="https://r-nimble.org/"> Nimble </a> ] .right-column9[  ] --- # Bayesian Logistic Regression using JAGS .left-column8[ - .hlb[Likelihood] $$ y_i \sim \text{Bernoulli}(\pi_i)\,, i= 1, \ldots, 400\,, $$ using logit link: $$ logit(\pi_i) = \beta_0 + \beta_1 age60_i + \beta_2 sbp140_i$$ - .hlb[Prior distributions] (weakly-informative): `$$\beta_0 \sim \mathcal{N}(0, 10^3),$$` `$$\beta_1 \sim \mathcal{N}(0, 10^3),$$` `$$\beta_2 \sim \mathcal{N}(0, 10^3).$$` ] .right-column9[ ``` r model_string <- " model { for (i in 1:N) { y[i] ~ dbern(pi[i]) logit(pi[i]) <- beta0 + beta1 * age60[i] + beta2 * sbp140[i] } # Priors for regression coefficients beta0 ~ dnorm(0, 0.001) beta1 ~ dnorm(0, 0.001) beta2 ~ dnorm(0, 0.001) }" ``` Check `S1-JAGS-heart_attack.Rmd` for the complete solution ] --- # Bayesian Logistic Regression using brms .left-column8[ - .hlb[Likelihood] $$ y_i \sim \text{Bernoulli}(\pi_i)\,, i= 1, \ldots, 400\,, $$ using logit link: $$ \text{logit}(\pi_i) = \beta_0 + \beta_1 \, \text{age60}_i + \beta_2 \, \text{sbp140}_i $$ - .hlb[Prior distributions] (weakly-informative): `$$\beta_0 \sim \mathcal{N}(0, 10^3),$$` `$$\beta_1 \sim \mathcal{N}(0, 10^3),$$` `$$\beta_2 \sim \mathcal{N}(0, 10^3).$$` ] .right-column9[ ``` r formula <- bf(y ~ age60 + sbp140, family = bernoulli()) # Fit the model using brms fit_brms <- brm(formula, data = data_hattack, prior = priors, chains = 3, # Number of MCMC chains iter = 5000, # Total number of iterations per chain warmup = 1000, # Number of iterations for warm-up thin = 1, # Thinning interval seed = 123, # Seed for reproducibility ) # Summary of the fitted model summary(fit) ``` Check `S1-brms-heart_attack.Rmd` for the complete solution ] --- # What we have learned so far - ### The biggest challenge for Bayesian inference has always been the .hlb[computational power] that more complex models require. -- - ### .hlb[MCMC methods allow computationally estimating posterior distributions] that are analytically intractable. -- - ### Today, there are many, many teams of researchers developing computational techniques for computing a posteriori distributions. --- # References .hlb[Books] - .hlb[Stan & R]: Lambert, B. (2018). A Student’s Guide to Bayesian Statistics. SAGE Publications. - .hlb[JAGS]: Plummer, M. (2019). JAGS User Manual Version 4.3.0. - .hlb[OpenBUGS]: Cowles, M. K. (2013). Applied Bayesian statistics: with R and OpenBUGS examples (Vol. 98). Springer Science & Business Media. - .hlb[WinBUGS]: Ntzoufras, I. (2011). Bayesian modeling using WinBUGS. John Wiley & Sons. - .hlb[Stan]: Andrew Gelman, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, Donald B. Rubin (2013). Bayesian Data Analysis. Chapman and Hall/CRC - .hlb[INLA]: Gómez-Rubio, V. (2020). Bayesian inference with INLA. CRC Press. <!-- - .hlb[Disease mapping]: Martínez-Beneito, M. A., & Botella-Rocamora, P. (2019). Disease mapping: from foundations to multidimensional modeling. CRC Press. --> <!-- - .hlb[Bayesian modelling]: Congdon, P. (2007). Bayesian statistical modelling. John Wiley & Sons. --> --- #References .hlb[Blogs] - http://wlm.userweb.mwn.de/R/wlmRcoda.htm - https://rpubs.com/FJRubio/IntroMCMC - https://darrenjw.wordpress.com/tag/mcmc/ - https://www.tweag.io/blog/2019-10-25-mcmc-intro1/ --- class: inverse, left, middle, animated, rotateInUpLeft </br> # ADIM: Bayesian GLM. MCMC methods ## Master's Degree in Data Analysis, Process Improvement and Decision Support Engineering </br> <font size="6"> Joaquín Martínez-Minaya, 2024-10-21 </font> </br> <a href="http://vabar.es/" style="color:white;" > VAlencia BAyesian Research Group </a> </br> <a href="https://smeg-bayes.org/ " style="color:white;"> Statistical Modeling Ecology Group </a> </br> <a href="https://giem.blogs.upv.es/" style="color:white;"> Grupo de Ingeniería Estadística Multivariante </a> </br> <a href="jmarmin@eio.upv.es" style="color:white;"> jmarmin@eio.upv.es </a>